In RodolfoPelinson/PredatorIsolationStochasticity: Data Analyses for Pelinson et al 2021B
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = TRUE)
These are the analyses of community size for the whole communities (results shown in the main papaer) and separetely to predatory and non-predatory insects.
These are the packages you will need to run this code:
library(lme4) # Version 1.1-23
library(emmeans) # Version 1.4.8
library(car) # Version 3.0-7
library(DHARMa) # Version 0.3.3.0
Whole Community
First, lets check what probability distribution should we choose using the most complex model we have:
mix_model_G <- lmer(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 +
(1|ID_SS2_SS3), REML = F)
mix_model_P <- glmer(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 +
(1|ID_SS2_SS3), family = "poisson")
mix_model_NB <- glmer.nb(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*
SS_SS2_SS3 + (1|ID_SS2_SS3))
resid_model_G <- simulateResiduals(mix_model_G)
plot(resid_model_G)
resid_model_P <- simulateResiduals(mix_model_P)
plot(resid_model_P)
resid_model_NB <- simulateResiduals(mix_model_NB)
plot(resid_model_NB)
AIC(mix_model_G,mix_model_P,mix_model_NB)
Looking at AIC values and the spread of residuals, it seems like negative binomial is the best option here.
So we can go further analysing the data:
mix_model_NB <- glmer.nb(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*
SS_SS2_SS3 + (1|ID_SS2_SS3))
round(Anova(mix_model_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_NB, list(pairwise ~ fish_SS2_SS3|SS_SS2_SS3),
adjust = "sidak")
emmeans(mix_model_NB, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3),
adjust = "sidak")
It seems that community size grows with time. But it grows larger in fishless ponds, and in higher isolation treatments.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(com_SS2_SS3_abundance~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(com_SS2_SS3_abundance)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_abundance[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_abundance~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Community Size", cex.lab = 1.3, line = 3)
title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75)
Ploting two surveys separetely
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS2)
#levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("transparent", "transparent"),3)
pch <- c(22,22,21,21,24,24) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Community Size", cex.lab = 1.3, line = 3)
title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75)
title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS3)
#levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("sienna3", "dodgerblue3"),3)
pch <- c(15,15,16,16,17,17) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Community Size", cex.lab = 1.3, line = 3)
title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75)
title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
Only Predatory Insects Community
First, lets load the necessary data:
data(com_SS2_SS3_predators_abundance)
Analysing the data:
mix_model_predators_NB <- glmer.nb(com_SS2_SS3_predators_abundance~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3),
control = glmerControl(optimizer = "bobyqa"))
round(Anova(mix_model_predators_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_predators_NB, list(pairwise ~ isolation_SS2_SS3), adjust = "sidak")
emmeans(mix_model_predators_NB, list(pairwise ~ fish_SS2_SS3|SS_SS2_SS3), adjust = "sidak")
It seems that when we only consider predatory insects, community size is affected by time (survey) presence of fish (which have a stronger effect in the third survey) and isolation.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(com_SS2_SS3_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(com_SS2_SS3_predators_abundance)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_predators_abundance[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2)
Ploting two surveys sseparetely
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS2)
#levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("transparent", "transparent"),3)
pch <- c(22,22,21,21,24,24) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2)
title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS3)
#levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("sienna3", "dodgerblue3"),3)
pch <- c(15,15,16,16,17,17) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2)
title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
Only Non-Predatory Insects (Herbivores and Detritivores) Community
First, lets load the necessary data:
data(com_SS2_SS3_non_predators_abundance)
Analysing the data:
mix_model_non_predators_NB <- glmer.nb(com_SS2_SS3_non_predators_abundance~fish_SS2_SS3*
isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3),
control = glmerControl(optimizer = "bobyqa"))
round(Anova(mix_model_non_predators_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_non_predators_NB, list(pairwise ~ SS_SS2_SS3|isolation_SS2_SS3),
adjust = "sidak")
It seems that we have similar patterns to when we consider the whole community.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1))
boxplot(com_SS2_SS3_non_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n")
mylevels <- levels(All)
#levelProportions <- summary(All)/length(com_SS2_SS3_non_predators_abundance)
col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(22,21,24,22,21,24,15,16,17,15,16,17)
for(i in 1:length(mylevels)){
x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_non_predators_abundance[All==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_non_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance", cex.lab = 1.3, line = 3)
title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75)
Ploting two surveys sseparetely
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS2)
#levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("transparent", "transparent"),3)
pch <- c(22,22,21,21,24,24) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance", cex.lab = 1.3, line = 3)
title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75)
title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1))
boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3,
outline = F, ylab = "", xlab = "",
at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n")
mylevels <- levels(fish_isolation_SS3)
#levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"])
col <- rep(c("sienna3", "dodgerblue3"),3)
bg <- rep(c("sienna3", "dodgerblue3"),3)
pch <- c(15,15,16,16,17,17) #o outro é 22,21,24
for(i in 1:length(mylevels)){
x<- c(1,5,2,6,3,7)[i]
thislevel <- mylevels[i]
thisvalues <- com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel]
# take the x-axis indices and add a jitter, proportional to the N in each level
myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3)
points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5)
}
boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3,
add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7),
lwd = 1.5, xaxt="n", yaxt="n")
axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1,
at =c(1,2,3,5,6,7), gap.axis = -10)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F )
axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE)
axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE)
title(ylab = "Abundance", cex.lab = 1.3, line = 3)
title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75)
title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
\
Herbivores and Detritivores VS Predators
Because we saw that community size increases with isolation, and we know from previous work that herbivores and detritivores are positively affected by isolation, we checked if herbivores and detritivores are indeed responsible for this increase.
First for the second Survey
#Second Survey
non_pred_barplot_SS2 <- tapply(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"], interaction(fish_SS2, isolation_SS2), sum)
non_pred_SS2 <- tapply(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"], interaction(isolation_SS2,fish_SS2), sum)
matrix_barplot_SS2 <- matrix(c(non_pred_SS2,non_pred_barplot_SS2), byrow = TRUE, ncol = 6, nrow = 2)
total_barplot_SS2 <- colSums(matrix_barplot_SS2)
percentage_barplot_SS2 <- matrix_barplot_SS2
percentage_barplot_SS2[1,] <- percentage_barplot_SS2[1,]/total_barplot_SS2
percentage_barplot_SS2[2,] <- percentage_barplot_SS2[2,]/total_barplot_SS2
barplot(percentage_barplot_SS2, space = c(0,0.5,0.5,1,0.5,0.5), xaxt="n", yaxt="n", col = c("lightgoldenrod2", "chartreuse4"))
axis(1,at = c(0.5,2,3.5,5.5,7,8.5), labels =rep(c("30 m","120 m","480 m"),2), tick = T, cex.axis = 0.9, gap.axis = 0)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.1, at =c(2,7), line = 1.5, tick = F )
axis(2, cex.axis = 1.2, gap.axis = 0, line = -0.5, tick = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1), labels = c("0%","20%","40%","60%","80%","100%"))
axis(2, cex.axis = 1.2, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1))
title(ylab = "Proportion of Abundance", cex.lab = 1.25, line = 3)
title(main = "Second Survey", cex.lab = 1.25, line = 2)
Third Survey
#Third Survey
non_pred_barplot_SS3 <- tapply(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"], interaction(fish_SS3, isolation_SS3), sum)
non_pred_SS3 <- tapply(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"], interaction(fish_SS3, isolation_SS3), sum)
matrix_barplot_SS3 <- matrix(c(non_pred_SS3,non_pred_barplot_SS3), byrow = TRUE, ncol = 6, nrow = 2)
total_barplot_SS3 <- colSums(matrix_barplot_SS3)
percentage_barplot_SS3 <- matrix_barplot_SS3
percentage_barplot_SS3[1,] <- percentage_barplot_SS3[1,]/total_barplot_SS3
percentage_barplot_SS3[2,] <- percentage_barplot_SS3[2,]/total_barplot_SS3
barplot(percentage_barplot_SS3, space = c(0,0.5,0.5,1,0.5,0.5), xaxt="n", yaxt="n", col = c("lightgoldenrod2", "chartreuse4"))
axis(1,at = c(0.5,2,3.5,5.5,7,8.5), labels =rep(c("30 m","120 m","480 m"),2), tick = T, cex.axis = 0.9, gap.axis = 0)
axis(1,labels = c("Fishless","Fish"), cex.axis = 1.1, at =c(2,7), line = 1.5, tick = F )
axis(2, cex.axis = 1.2, gap.axis = 0, line = -0.5, tick = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1), labels = c("0%","20%","40%","60%","80%","100%"))
axis(2, cex.axis = 1.2, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1))
title(ylab = "Proportion of Abundance", cex.lab = 1.25, line = 3)
title(main = "Third Survey", cex.lab = 1.25, line = 2)
It seems that herbivores and detritivores are indeed responsible for the patterns observed in community size. More importantly, this effect is stronger in more isolated ponds.
RodolfoPelinson/PredatorIsolationStochasticity documentation built on Feb. 21, 2022, 12:03 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(echo = TRUE)
These are the analyses of community size for the whole communities (results shown in the main papaer) and separetely to predatory and non-predatory insects.
These are the packages you will need to run this code:
library(lme4) # Version 1.1-23 library(emmeans) # Version 1.4.8 library(car) # Version 3.0-7 library(DHARMa) # Version 0.3.3.0
Whole Community
First, lets check what probability distribution should we choose using the most complex model we have:
mix_model_G <- lmer(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), REML = F) mix_model_P <- glmer(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), family = "poisson") mix_model_NB <- glmer.nb(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3* SS_SS2_SS3 + (1|ID_SS2_SS3)) resid_model_G <- simulateResiduals(mix_model_G) plot(resid_model_G) resid_model_P <- simulateResiduals(mix_model_P) plot(resid_model_P) resid_model_NB <- simulateResiduals(mix_model_NB) plot(resid_model_NB) AIC(mix_model_G,mix_model_P,mix_model_NB)
Looking at AIC values and the spread of residuals, it seems like negative binomial is the best option here.
So we can go further analysing the data:
mix_model_NB <- glmer.nb(com_SS2_SS3_abundance~fish_SS2_SS3*isolation_SS2_SS3* SS_SS2_SS3 + (1|ID_SS2_SS3)) round(Anova(mix_model_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_NB, list(pairwise ~ fish_SS2_SS3|SS_SS2_SS3), adjust = "sidak") emmeans(mix_model_NB, list(pairwise ~ isolation_SS2_SS3|SS_SS2_SS3), adjust = "sidak")
It seems that community size grows with time. But it grows larger in fishless ponds, and in higher isolation treatments.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(com_SS2_SS3_abundance~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(com_SS2_SS3_abundance) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_abundance[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_abundance~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Community Size", cex.lab = 1.3, line = 3) title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75)
Ploting two surveys separetely
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS2) #levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("transparent", "transparent"),3) pch <- c(22,22,21,21,24,24) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Community Size", cex.lab = 1.3, line = 3) title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75) title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS3) #levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("sienna3", "dodgerblue3"),3) pch <- c(15,15,16,16,17,17) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Community Size", cex.lab = 1.3, line = 3) title(ylab = "(Total Abundance)", cex.lab = 1.3, line = 1.75) title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
Only Predatory Insects Community
First, lets load the necessary data:
data(com_SS2_SS3_predators_abundance)
Analysing the data:
mix_model_predators_NB <- glmer.nb(com_SS2_SS3_predators_abundance~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = glmerControl(optimizer = "bobyqa")) round(Anova(mix_model_predators_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_predators_NB, list(pairwise ~ isolation_SS2_SS3), adjust = "sidak") emmeans(mix_model_predators_NB, list(pairwise ~ fish_SS2_SS3|SS_SS2_SS3), adjust = "sidak")
It seems that when we only consider predatory insects, community size is affected by time (survey) presence of fish (which have a stronger effect in the third survey) and isolation.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(com_SS2_SS3_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(com_SS2_SS3_predators_abundance) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_predators_abundance[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2)
Ploting two surveys sseparetely
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS2) #levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("transparent", "transparent"),3) pch <- c(22,22,21,21,24,24) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2) title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,250), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS3) #levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("sienna3", "dodgerblue3"),3) pch <- c(15,15,16,16,17,17) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance (Predators)", cex.lab = 1.3, line = 2) title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
Only Non-Predatory Insects (Herbivores and Detritivores) Community
First, lets load the necessary data:
data(com_SS2_SS3_non_predators_abundance)
Analysing the data:
mix_model_non_predators_NB <- glmer.nb(com_SS2_SS3_non_predators_abundance~fish_SS2_SS3* isolation_SS2_SS3*SS_SS2_SS3 + (1|ID_SS2_SS3), control = glmerControl(optimizer = "bobyqa")) round(Anova(mix_model_non_predators_NB, test.statistic = "Chisq"),3)
Now pairwise differences:
emmeans(mix_model_non_predators_NB, list(pairwise ~ SS_SS2_SS3|isolation_SS2_SS3), adjust = "sidak")
It seems that we have similar patterns to when we consider the whole community.
Lets plot it:
Ploting two surveys together
par(cex = 0.75, mar = c(4,4,0.1,0.1)) boxplot(com_SS2_SS3_non_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", gap.axis = -100, yaxt="n") mylevels <- levels(All) #levelProportions <- summary(All)/length(com_SS2_SS3_non_predators_abundance) col <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3)) bg <- c(rep("transparent",3), rep("transparent",3),rep("sienna3",3), rep("dodgerblue3",3)) pch <- c(22,21,24,22,21,24,15,16,17,15,16,17) for(i in 1:length(mylevels)){ x<- c(1,2,3,5,6,7,1,2,3,5,6,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_non_predators_abundance[All==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_non_predators_abundance~isolation_SS2_SS3*fish_SS2_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance", cex.lab = 1.3, line = 3) title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75)
Ploting two surveys sseparetely
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS2) #levelProportions <- summary(fish_isolation_SS2)/length(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("transparent", "transparent"),3) pch <- c(22,22,21,21,24,24) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"][fish_isolation_SS2==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]~isolation_SS2*fish_SS2, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance", cex.lab = 1.3, line = 3) title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75) title(main = "Second Survey", cex.lab = 1.3, line = 0.5)
par(cex = 0.75, mar = c(4,4,1.5,0.1)) boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3, outline = F, ylab = "", xlab = "", at = c(1,2,3,5,6,7),ylim = c(0,1400), lwd = 1.5, col = "transparent", xaxt="n", yaxt="n") mylevels <- levels(fish_isolation_SS3) #levelProportions <- summary(fish_isolation_SS3)/length(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"]) col <- rep(c("sienna3", "dodgerblue3"),3) bg <- rep(c("sienna3", "dodgerblue3"),3) pch <- c(15,15,16,16,17,17) #o outro é 22,21,24 for(i in 1:length(mylevels)){ x<- c(1,5,2,6,3,7)[i] thislevel <- mylevels[i] thisvalues <- com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"][fish_isolation_SS3==thislevel] # take the x-axis indices and add a jitter, proportional to the N in each level myjitter <- jitter(rep(x, length(thisvalues)), amount=0.3) points(myjitter, thisvalues, pch=pch[i], col=col[i], bg = bg[i] , cex = 1.5, lwd = 1.5) } boxplot(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"]~isolation_SS3*fish_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n", yaxt="n") axis(1,labels = c("30","120", "480","30","120", "480"), cex.axis = 1.1, at =c(1,2,3,5,6,7), gap.axis = -10) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.3, at =c(2,6), line = 1.5, tick = F ) axis(2, cex.axis = 0.8, gap.axis = 0, line = -0.5, tick = FALSE) axis(2, cex.axis = 0.8, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE) title(ylab = "Abundance", cex.lab = 1.3, line = 3) title(ylab = "(Herbivores and Detritivores)", cex.lab = 1.3, line = 1.75) title(main = "Third Survey", cex.lab = 1.3, line = 0.5)
\
Herbivores and Detritivores VS Predators
Because we saw that community size increases with isolation, and we know from previous work that herbivores and detritivores are positively affected by isolation, we checked if herbivores and detritivores are indeed responsible for this increase.
First for the second Survey
#Second Survey non_pred_barplot_SS2 <- tapply(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "2"], interaction(fish_SS2, isolation_SS2), sum) non_pred_SS2 <- tapply(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "2"], interaction(isolation_SS2,fish_SS2), sum) matrix_barplot_SS2 <- matrix(c(non_pred_SS2,non_pred_barplot_SS2), byrow = TRUE, ncol = 6, nrow = 2) total_barplot_SS2 <- colSums(matrix_barplot_SS2) percentage_barplot_SS2 <- matrix_barplot_SS2 percentage_barplot_SS2[1,] <- percentage_barplot_SS2[1,]/total_barplot_SS2 percentage_barplot_SS2[2,] <- percentage_barplot_SS2[2,]/total_barplot_SS2 barplot(percentage_barplot_SS2, space = c(0,0.5,0.5,1,0.5,0.5), xaxt="n", yaxt="n", col = c("lightgoldenrod2", "chartreuse4")) axis(1,at = c(0.5,2,3.5,5.5,7,8.5), labels =rep(c("30 m","120 m","480 m"),2), tick = T, cex.axis = 0.9, gap.axis = 0) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.1, at =c(2,7), line = 1.5, tick = F ) axis(2, cex.axis = 1.2, gap.axis = 0, line = -0.5, tick = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1), labels = c("0%","20%","40%","60%","80%","100%")) axis(2, cex.axis = 1.2, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1)) title(ylab = "Proportion of Abundance", cex.lab = 1.25, line = 3) title(main = "Second Survey", cex.lab = 1.25, line = 2)
Third Survey
#Third Survey non_pred_barplot_SS3 <- tapply(com_SS2_SS3_non_predators_abundance[SS_SS2_SS3 == "3"], interaction(fish_SS3, isolation_SS3), sum) non_pred_SS3 <- tapply(com_SS2_SS3_predators_abundance[SS_SS2_SS3 == "3"], interaction(fish_SS3, isolation_SS3), sum) matrix_barplot_SS3 <- matrix(c(non_pred_SS3,non_pred_barplot_SS3), byrow = TRUE, ncol = 6, nrow = 2) total_barplot_SS3 <- colSums(matrix_barplot_SS3) percentage_barplot_SS3 <- matrix_barplot_SS3 percentage_barplot_SS3[1,] <- percentage_barplot_SS3[1,]/total_barplot_SS3 percentage_barplot_SS3[2,] <- percentage_barplot_SS3[2,]/total_barplot_SS3 barplot(percentage_barplot_SS3, space = c(0,0.5,0.5,1,0.5,0.5), xaxt="n", yaxt="n", col = c("lightgoldenrod2", "chartreuse4")) axis(1,at = c(0.5,2,3.5,5.5,7,8.5), labels =rep(c("30 m","120 m","480 m"),2), tick = T, cex.axis = 0.9, gap.axis = 0) axis(1,labels = c("Fishless","Fish"), cex.axis = 1.1, at =c(2,7), line = 1.5, tick = F ) axis(2, cex.axis = 1.2, gap.axis = 0, line = -0.5, tick = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1), labels = c("0%","20%","40%","60%","80%","100%")) axis(2, cex.axis = 1.2, gap.axis = 0, line = 0, tick = TRUE, labels = FALSE, at = c(0.0,0.2,0.4,0.6,0.8,1)) title(ylab = "Proportion of Abundance", cex.lab = 1.25, line = 3) title(main = "Third Survey", cex.lab = 1.25, line = 2)
It seems that herbivores and detritivores are indeed responsible for the patterns observed in community size. More importantly, this effect is stronger in more isolated ponds.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.